Qu et al V2

1 Functions

2 Data

# read processed data:
SCC1_cancer  = readRDS(file = "./Data/cervical_cancer_data/Qiu et al/SCC1/SCC1_cancer.RDS")
SCC2_cancer  = readRDS(file = "./Data/cervical_cancer_data/Qiu et al/SCC2/SCC2_cancer_processed.RDS")
SCC3_cancer = readRDS(file = "./Data/cervical_cancer_data/Qiu et al/SCC3/SCC3_cancer_processed.RDS")
SCC4_cancer = readRDS(file = "./Data/cervical_cancer_data/Qiu et al/SCC4/SCC4_cancer_processed.RDS")
SCC5_cancer  = readRDS(file = "./Data/cervical_cancer_data/Qiu et al/SCC5/SCC5_cancer_processed.RDS")

scc.big <- merge(SCC1_cancer, y = c(SCC2_cancer,SCC3_cancer, SCC4_cancer,SCC5_cancer), add.cell.ids = c("SCC1","SCC2", "SCC3", "SCC4","SCC5"), project = "SCC")
scc.big$orig.ident =  sapply(X = strsplit(colnames(scc.big), split = "_"), FUN = "[", 1)


VlnPlot(object = scc.big,features = "MYB",group.by = "orig.ident",slot = "data", assay = "RNA")

NA
NA

myb_data = FetchData(object = scc.big,vars = c("MYB","orig.ident"),slot = "data", assay = "RNA")
myb_data %>%   group_by(orig.ident) %>% 
  summarise(counts = sum(MYB > 0, na.rm = TRUE))

3 Ppatients filter

scc_myb_patients<- subset(scc.big, subset = orig.ident %in% c("SCC3","SCC4","SCC5"))

pdf(file = "./Figures/SCC_MYB_HPV_boxplot_all_patients.pdf",height = 4)
p
dev.off()

null device 
          1 

4 Boxplot by patient

scc_myb_patients<- subset(scc.big, subset = orig.ident %in% c("SCC2","SCC3","SCC4","SCC5"))

df = FetchData(object = scc_myb_patients,vars = c("orig.ident","MYB","hpv_positive")) %>% group_by(orig.ident,hpv_positive) %>% summarise(
  myb = mean(MYB),.groups = "drop_last")
stat.test <- df %>% group_by("hpv_positive") %>% 
  wilcox_test(myb ~ hpv_positive)


df = reshape2::dcast(df, orig.ident ~ hpv_positive)

Using myb as value column: use value.var to override.

p = ggpaired(df, cond1  = "positive",cond2   = "negative",palette = "jco",
            add = "jitter", line.color = "gray", line.size = 0.4,id = "orig.ident", fill = "condition") + theme_minimal()+
  # stat_compare_means(method = "wilcox.test",comparisons = list(c("positive","negative")),paired = F)+
  stat_summary(fun.data = function(x) data.frame(y=max(x)*1.2, label = paste("Mean=",round(mean(x),digits = 2))), geom="text")+ylab("MYB")+ stat_pvalue_manual(stat.test,y.position = 0.2, label = "Wilcox, p = {p}",remove.bracket = F)

Warning: `gather_()` was deprecated in tidyr 1.2.0.
Please use `gather()` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

p

pdf(file = "./Figures/SCC_HPV_MYB_boxplot_bySample.pdf",width = 8)
p
dev.off()

null device 
          1 

5 HPV DEG boxplots

genes = c("ANLN", "TUBB6", "AXL", "IFT122", "GFM1", "MYB", "JAG1")
myb_vs_hpv = FetchData(object = scc_myb_patients,vars = c("hpv_positive",genes))
df = reshape2::melt(myb_vs_hpv,value.name = "logTPM") %>% dplyr::rename(gene = variable)

Using hpv_positive as id variables

stat.test <- df %>%
    group_by(gene) %>%
  wilcox_test(logTPM ~ hpv_positive) %>%
  mutate(y.position = 5)

stat.test

stat.test <- stat.test %>% 
  add_xy_position(x = "gene", dodge = 0.8)

p = ggboxplot(
  df,
  x = "gene",
  y = "logTPM",
  color = "hpv_positive",
  palette = "jco",
  add = "jitter"
)+ stat_pvalue_manual(stat.test,y.position = 4, label = "p = {p}",remove.bracket = T)

p

pdf(file = "./Figures/SCC_HPV_genes_all_patients.pdf",width = 8)
p
dev.off()

metagenes_violin_compare.2 = function(dataset,prefix = "",pre_on = c("OSI","NT"),axis.text.x = 11,test = "t.test", programs = c("Hypoxia","TNFa","Cell_cycle"),return_list = F,combine_patients = F){
  require(facefuns)
  plt.lst = list()
  if(combine_patients){
    genes_by_tp = FetchData(object = dataset,vars =  c("treatment",programs)) %>% filter(treatment %in% pre_on)  %>% as.data.frame() #mean expression
    formula <- as.formula( paste("c(", paste(programs, collapse = ","), ")~ treatment ") )
    
    #plot and split by patient:   
    stat.test = compare_means(formula = formula ,data = genes_by_tp,method = test,p.adjust.method = "fdr")%>% # Add pairwise comparisons p-value
      dplyr::filter(group1 == pre_on[1] & group2 == pre_on[2])  #filter for pre vs on treatment only
    
    stat.test$p.format =stat.test$p.adj #modift 0 pvalue to be lowest possible float
    stat.test$p.format[!stat.test$p.format == 0 ] <- paste("=",stat.test$p.format[!stat.test$p.format == 0 ])
    stat.test$p.format[stat.test$p.format == 0 ] <- paste("<",.Machine$double.xmin %>% signif(digits = 3))
    
    
    genes_by_tp = reshape2::melt(genes_by_tp, id.vars = c("treatment"),value.name = "score")
    plt = ggplot(genes_by_tp, aes(x = variable, y = score,fill = treatment)) + geom_split_violin(scale = 'width')+ 
      geom_boxplot(width = 0.25, notch = FALSE, notchwidth = .4, outlier.shape = NA, coef=0)+
      ylim(min(genes_by_tp$score),max(genes_by_tp$score)*1.25)
    plt = plt +stat_pvalue_manual(stat.test, label = "asd = {p.signif}", label.size = 7,  #add p value
                                  y.position = 1,inherit.aes = F,size = 3.3,x = ".y.") # set position at the top value
    return(plt)
  }
  
  
  
  
  
  if (return_list) {
    return(plt.lst)
  }
}

scc_myb_patients$treatment = scc_myb_patients$hpv_positive %>% gsub(pattern = "negative",replacement = "HPV-negative")%>% gsub(pattern = "positive",replacement = "HPV-positive")
scc_myb_patients@meta.data[["treatment"]] = factor(scc_myb_patients$treatment, levels = c("HPV-positive", "HPV-negative"))

p = metagenes_violin_compare.2(dataset = scc_myb_patients, prefix = "patient",pre_on = c("HPV-negative","HPV-positive"),test = "wilcox.test",programs = genes, return_list = F,combine_patients = T) +scale_y_continuous(limits = c(0,1.5)) + labs(fill = "")+ylab("LogTPM")+xlab("Gene")+theme(axis.title=element_text(size=14))+ scale_fill_manual(values=c("#F8766D", "cyan3"))

Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

p

Warning: Removed 398 rows containing non-finite values (stat_ydensity).
Warning: Removed 398 rows containing non-finite values (stat_boxplot).

pdf(file = "./Figures/SCC_HPV_Top_violin.pdf",width = 10,height = 5)
p

Warning: Removed 392 rows containing non-finite values (stat_ydensity).
Warning: Removed 392 rows containing non-finite values (stat_boxplot).

dev.off()

null device 
          1 

myplots <- list()  # new empty list

for (patient_name in c("SCC2", "SCC3", "SCC4", "SCC5")) {
  patient_data = subset(x = scc_myb_patients, subset = orig.ident == patient_name)
  df  = FetchData(object = patient_data,
                  vars = c("hpv_positive", "MYB_positive")) %>% droplevels()
  test = fisher_test(table(df))
  p = ggbarstats(
    df,
    MYB_positive,
    hpv_positive,
    results.subtitle = FALSE,
    subtitle = paste0("Fisher's exact test", ", p-value = ",
                      test$p)
  ,title = patient_name)
  myplots[[patient_name]] <- p  # add each plot into plot list

}
ggarrange(plotlist = myplots,ncol = 4,nrow = 1)

deg <-
  FindMarkers(
    scc_myb_patients,
    ident.1 = "positive",
    ident.2 = "negative",
    features = features,
    densify = T,
    assay = "RNA",
    logfc.threshold = 0.1,
    min.pct = 0.1,
    only.pos = F,
    mean.fxn = function(x) {
      return(log(rowMeans(x) + 1, base = 2)) # change func to calculate logFC in log space data (default to exponent data)
    }
  )

  |                                                  | 0 % ~calculating  
  |+                                                 | 1 % ~14s          
  |++                                                | 2 % ~13s          
  |++                                                | 3 % ~14s          
  |+++                                               | 4 % ~13s          
  |+++                                               | 5 % ~13s          
  |++++                                              | 6 % ~13s          
  |++++                                              | 7 % ~13s          
  |+++++                                             | 8 % ~13s          
  |+++++                                             | 9 % ~13s          
  |++++++                                            | 10% ~13s          
  |++++++                                            | 11% ~13s          
  |+++++++                                           | 12% ~13s          
  |+++++++                                           | 13% ~13s          
  |++++++++                                          | 14% ~13s          
  |++++++++                                          | 15% ~12s          
  |+++++++++                                         | 16% ~12s          
  |+++++++++                                         | 18% ~12s          
  |++++++++++                                        | 19% ~12s          
  |++++++++++                                        | 20% ~12s          
  |+++++++++++                                       | 21% ~12s          
  |+++++++++++                                       | 22% ~12s          
  |++++++++++++                                      | 23% ~12s          
  |++++++++++++                                      | 24% ~11s          
  |+++++++++++++                                     | 25% ~12s          
  |+++++++++++++                                     | 26% ~11s          
  |++++++++++++++                                    | 27% ~11s          
  |++++++++++++++                                    | 28% ~11s          
  |+++++++++++++++                                   | 29% ~11s          
  |+++++++++++++++                                   | 30% ~11s          
  |++++++++++++++++                                  | 31% ~11s          
  |++++++++++++++++                                  | 32% ~10s          
  |+++++++++++++++++                                 | 33% ~10s          
  |++++++++++++++++++                                | 34% ~10s          
  |++++++++++++++++++                                | 35% ~10s          
  |+++++++++++++++++++                               | 36% ~10s          
  |+++++++++++++++++++                               | 37% ~10s          
  |++++++++++++++++++++                              | 38% ~09s          
  |++++++++++++++++++++                              | 39% ~09s          
  |+++++++++++++++++++++                             | 40% ~09s          
  |+++++++++++++++++++++                             | 41% ~09s          
  |++++++++++++++++++++++                            | 42% ~09s          
  |++++++++++++++++++++++                            | 43% ~09s          
  |+++++++++++++++++++++++                           | 44% ~08s          
  |+++++++++++++++++++++++                           | 45% ~08s          
  |++++++++++++++++++++++++                          | 46% ~08s          
  |++++++++++++++++++++++++                          | 47% ~08s          
  |+++++++++++++++++++++++++                         | 48% ~08s          
  |+++++++++++++++++++++++++                         | 49% ~08s          
  |++++++++++++++++++++++++++                        | 51% ~08s          
  |++++++++++++++++++++++++++                        | 52% ~07s          
  |+++++++++++++++++++++++++++                       | 53% ~07s          
  |+++++++++++++++++++++++++++                       | 54% ~07s          
  |++++++++++++++++++++++++++++                      | 55% ~07s          
  |++++++++++++++++++++++++++++                      | 56% ~07s          
  |+++++++++++++++++++++++++++++                     | 57% ~07s          
  |+++++++++++++++++++++++++++++                     | 58% ~06s          
  |++++++++++++++++++++++++++++++                    | 59% ~06s          
  |++++++++++++++++++++++++++++++                    | 60% ~06s          
  |+++++++++++++++++++++++++++++++                   | 61% ~06s          
  |+++++++++++++++++++++++++++++++                   | 62% ~06s          
  |++++++++++++++++++++++++++++++++                  | 63% ~06s          
  |++++++++++++++++++++++++++++++++                  | 64% ~06s          
  |+++++++++++++++++++++++++++++++++                 | 65% ~05s          
  |+++++++++++++++++++++++++++++++++                 | 66% ~05s          
  |++++++++++++++++++++++++++++++++++                | 67% ~05s          
  |+++++++++++++++++++++++++++++++++++               | 68% ~05s          
  |+++++++++++++++++++++++++++++++++++               | 69% ~05s          
  |++++++++++++++++++++++++++++++++++++              | 70% ~05s          
  |++++++++++++++++++++++++++++++++++++              | 71% ~05s          
  |+++++++++++++++++++++++++++++++++++++             | 72% ~05s          
  |+++++++++++++++++++++++++++++++++++++             | 73% ~04s          
  |++++++++++++++++++++++++++++++++++++++            | 74% ~04s          
  |++++++++++++++++++++++++++++++++++++++            | 75% ~04s          
  |+++++++++++++++++++++++++++++++++++++++           | 76% ~04s          
  |+++++++++++++++++++++++++++++++++++++++           | 77% ~04s          
  |++++++++++++++++++++++++++++++++++++++++          | 78% ~04s          
  |++++++++++++++++++++++++++++++++++++++++          | 79% ~03s          
  |+++++++++++++++++++++++++++++++++++++++++         | 80% ~03s          
  |+++++++++++++++++++++++++++++++++++++++++         | 81% ~03s          
  |++++++++++++++++++++++++++++++++++++++++++        | 82% ~03s          
  |++++++++++++++++++++++++++++++++++++++++++        | 84% ~03s          
  |+++++++++++++++++++++++++++++++++++++++++++       | 85% ~03s          
  |+++++++++++++++++++++++++++++++++++++++++++       | 86% ~02s          
  |++++++++++++++++++++++++++++++++++++++++++++      | 87% ~02s          
  |++++++++++++++++++++++++++++++++++++++++++++      | 88% ~02s          
  |+++++++++++++++++++++++++++++++++++++++++++++     | 89% ~02s          
  |+++++++++++++++++++++++++++++++++++++++++++++     | 90% ~02s          
  |++++++++++++++++++++++++++++++++++++++++++++++    | 91% ~02s          
  |++++++++++++++++++++++++++++++++++++++++++++++    | 92% ~01s          
  |+++++++++++++++++++++++++++++++++++++++++++++++   | 93% ~01s          
  |+++++++++++++++++++++++++++++++++++++++++++++++   | 94% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 95% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 96% ~01s          
  |+++++++++++++++++++++++++++++++++++++++++++++++++ | 97% ~01s          
  |+++++++++++++++++++++++++++++++++++++++++++++++++ | 98% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 99% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=16s  

intersect(
  acc_deg %>% filter(.$fdr<0.05) %>% filter(.$avg_log2FC>0)   %>% arrange(dplyr::desc(.$avg_log2FC))  %>% head(400) %>% rownames(),
  deg %>% filter(.$fdr<0.05) %>% filter(.$avg_log2FC>0) %>% arrange(dplyr::desc(.$avg_log2FC)) %>% head(400) %>% rownames()
) %>% cat(sep = "\n")

LSM5
ETS2
TMPO
PSME2
TMX1
MSH6
PHB
SNRPD1
MCM3
RNF168
MCM5
MT1E
FBLN1
TXN
TUBA1B
MT2A
STMN1

# HNSC OS p=0.022
intersect(
   acc_deg %>% filter(.$fdr<0.05) %>% arrange(dplyr::desc(.$avg_log2FC))%>% filter(.$avg_log2FC>0) %>% head(200) %>% rownames(),
   deg %>% filter(.$fdr<0.05) %>% arrange(dplyr::desc(.$avg_log2FC))%>% filter(.$avg_log2FC>0) %>% head(200) %>% rownames()
 ) %>% cat(sep = "\n")

TMPO
PSME2
TMX1
SNRPD1
MCM3
MT1E
MT2A
STMN1

intersect(
   acc_deg %>% filter(.$fdr<0.05) %>% arrange(dplyr::desc(.$fdr))%>% filter(.$avg_log2FC>0) %>% head(200) %>% rownames(),
   deg %>% filter(.$fdr<0.05) %>% arrange(dplyr::desc(.$fdr))%>% filter(.$avg_log2FC>0) %>% head(200) %>% rownames()
 ) %>% cat(sep = "\n")

LGALS1
GAPDH
TUBA1C
CAPG
TAGLN2

intersect( 
   acc_deg %>% filter(.$fdr<0.05) %>% arrange(dplyr::desc(.$fdr))%>% filter(.$avg_log2FC>0) %>% head(400) %>% rownames(),
   deg %>% filter(.$fdr<0.05) %>% arrange(dplyr::desc(.$fdr))%>% filter(.$avg_log2FC>0) %>% head(400) %>% rownames()
 ) %>% cat(sep = "\n")

LGALS1
GAPDH
NDUFS7
MTHFD2
FDFT1
ETS2
KLF10
TUBA1C
CAPG
TAGLN2

# CECS RFS p=0.014
deg %>% filter(.$fdr<0.05) %>% filter(.$avg_log2FC>0) %>% arrange(dplyr::desc(.$avg_log2FC)) %>% rownames() %>% head(20) %>% cat(sep = "\n")

KRT5
IGFBP2
AQP3
HES1
ATP1B3
MT1X
NASP
DEK
SCPEP1
ACTL6A
FABP5
KRT15
GOLIM4
STMN1
ID1
DSC3
PRKDC
LMO4
SELENOP
NUCKS1

acc_deg[deg %>% filter(.$fdr<0.05) %>% filter(.$avg_log2FC>0) %>% arrange(dplyr::desc(.$avg_log2FC)) %>% rownames() %>% head(20),]